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FIRST QUARTER FIRST YEAR PROGRESS REPORT 

This report covers technical progress during the first quarter of the first year of NASA 
Sun-Earth Connections Theory Program (SECTP) contract “The Structure and Dynamics of 
the Solar Corona and Inner Heliosphere,” NA'S5-99188/between NASA and Science 
Applications International Corporation, and covers dtej^erida August 16, 1999 to November 15, 
1999. Under this contract SAIC and the University of California, Irvine (UCI) have conducted 
research into theoretical modeling of active regions, the solar corona, and the inner heliosphere, 
using the MHD model. ■ ~- 

In the following sections we summarize our progress during this reporting period. Full 
descriptions of our work can be found in the cited publications, a few of which are attached to 
this report. 

Publication on the August 11, 1999 Eclipse Prediction 

Our publication for the proceedings of the conference “The Last Total Solar Eclipse of 
the Millennium,” which was held in Istanbul, Turkey, August 13-15, 1999, has been completed 
(Mikic et al. 1999), and is included in the Appendix of this report. It discusses our 3D MHD 
modeling of the total solar eclipse that occurred on August 11, 1999, and a comparison between 
our model and observations. 

Publication on the Modeling of the Solar Wind and Transition Region 

Roberto Lionello presented a paper on “Magnetohydrodynamics of the Solar Corona and 
the Transition Region,” at the Meeting on Magnetic Fields and Solar Processes, which was 
held in Florence, Italy, from September 12-18, 1999. The paper has appeared in the 
proceedings of this conference (Lionello, Linker, & Mikic 1999), and is included in the 
Appendix. An extended paper on this work is in preparation (Lionello, Linker, & Mikic 2000). 

A New 3D MHD Code on Unstructured Tetrahedral Grids 

We have developed a new algorithm for the solution of the time-dependent, resistive 
MHD equations in three-dimensional spatial domains of arbitrary shape and connectivity. The 
algorithm uses a finite-volume approach based on a primary grid of tetrahedral cells, and a 
secondary dual median grid. The resulting discrete operators preserve the annihilation 
properties of the divergence and curl, and lead to a compact, self-adjoint formulation. The 
discrete operators also minimize the same functionals as the original differential operators. The 
semi-implicit method of time advancement is used. The capability of mesh refinement and 
coarsening has also been, built (but has not yet been tested). The algorithm is being 

C 

implemented in F90 and MPI, and can be used on either single processor or massively parallel 
computers. 

For initial testing purposes, we have used MH4D to solve the magnetic resistive equation 
and the scalar and vector advection equations. The flexibility and parallelism of MH4D is 
illustrated in Figure 1. We consider three of the many possible domain configurations (a cube, 
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Figure 1. Decomposition of various spatial grids among 5 processors for the MH4D code. The code solves 
the MHD equations on an unstructured tetrahedral grid of arbitrary shape and connectivity. We show how 
the mesh is split between 5 processors for three illustrative domain configurations (a cube, a cylinder, and a 
sphere). 
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a cylinder, and a sphere), showing how the tetrahedral mesh is split between 5 processors. This 
code is presently being extended to solve the full MHD equations, and will be used in the future 
on a variety of solar and other MHD problems. The presentation “MH4D: A New Algorithm 
for Three-Dimensional MHD on an Unstructured Tetrahedral Grid,” by R. Lionello and D. D. 
Schnack, was presented at the American Physical Society/Division of Plasma Physics Meeting, 
held in Seattle, WA, Nov. 15-19, 1999. 

Comparing MHD Simulations of the Inner Heliosphere with Spacecraft Magnetic 
Field Measurements 

Using the observed line-of-sight component of the photospheric magnetic field as input to 
a 3-dimensional MHD model, we have simulated the corona and heliosphere (out to 5 AU) for 
the time period known as “Whole Sun Month,” August-September, 1996, which occurred 
shortly after the minimum of solar cycle 22. By flying simulated trajectories through the model, 
we can compare directly these results with interplanetary spacecraft observations. In Figure 2 
we make such a comparison, using observations taken by the WIND and Ulysses spacecraft, 
which were widely separated in heliocentric distance, latitude, and longitude. Prior to this time, 
Ulysses was desc'ending from higher to lower latitudes, and this rotation was the first time that 
Ulysses encountered significantly variable wind. Ulysses was located on the opposite side of 
the Sun from Earth, at a heliocentric distance of 4.25 AU and heliographic latitude of - N 28 
degrees. In contrast, WIND was located in the ecliptic plane at 1 AU. The left column 
compares speed, scaled number density, and temperature at Ulysses, while the right column 
make the same comparison at WIND. In each panel, observations are colored blue and 
simulation results are red. At Ulysses, where the underlying profile consisted of a simple, 
single-stream pattern, the model appears to reproduce the essential features of the large-scale 
variations, i.e., the interaction region and the expansion wave. However, the interaction region is 
not as steep as the observations indicate, and the model does not mimic fluctuations on scales 
less than a few days. Nevertheless, when we consider that the measurements of the line-of-sight 
photospheric magnetic field used to drive the solution were always separated by ~ 180 degrees 

from the actual solar wind that Ulysses observed, the close agreement suggests that relatively 
little evolution of the corona took place, at least at latitudes of ~ 30 degrees. At WIND, there is 
a more complex pattern of variations that are only partially modeled by the MHD solution. The 
model predicts two fast (> 500 km/s) streams. Of these, one is found in the observations (at 
day 255), and the other (at day 240) might be related to the fast stream observed on day 243, 
although it is considerably broader. There are other minor perturbations in the model; however, 
none of them can be reliably matched with WIND observations. Overall, however, this 
comparison indicates that the model has reproduced the essential large-scale features of the 
heliosphere during this time period. 
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Figure 2. Comparison of Ulysses and WIND spacecraft measurements (blue) with 3D MHD simulations 
(red) during the Whole Sun Month time interval (August-September 1996). The Ulysses comparison is 
shown in the left panels; the WIND comparison is shown in the right panels. 
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Progress in Parallelizing the MAS Code 

Optimization of Inter-Process Communication 

The entire code has been carefully tested in order to avoid superfluous MPI 
communication calls. The size of the data sets to be sent/received in each communication call 
has been kept as large as possible by lumping together different arrays to maximize the 
communication efficiency. 

In order to improve further the communication efficiency, different MPI communication 
protocols (that can be selected at run time via an input' switch) have been introduced. In 
particular, the code can use “persistent communication routines” to reduce the overhead of 
repeated communication calls. Preliminary testing indicates a savings of up to 20% with respect 
the a more conventional coupled send/receive scheme. Since this performance-improvement can 
be problem dependent, further testing with different mesh sizes is required. 

To better assess code performance, a measure of the latency time (time “wasted” while 
waiting for the processors to become synchronized) has been introduced. In the tests 
conducted so far, the latency time appears to be at least an order of magnitude smaller than the 
communication time. 

Code Benchmarking 

Careful benchmarking of the parallel version vs. the serial MAS code has been conducted 
for different test problems. The parallel version, when run on a single processor, produces 
restart files that agree perfectly (with no difference error) with the serial version. When run 
with more than one processor, the agreement with the serial version is maintained within 
machine roundoff accuracy. The origin of this difference has been tracked down to the non- 
associative nature of floating point additions. 

Input-Controlled Problem Size 

The scalar MAS code and the first parallel versions that have been derived from it do not 
allow the user to change the problem size (i.e., mesh dimensions) at run time. To change the 
problem size requires a re-compilation of the code. A new version has been developed that 
allows the choice of the grid size from the input file. As an important consequence, this 
upgrade removes the constraint that required the same number of mesh point on each processor. 

Full testing of this new feature will be performed after this version is integrated with new I/O 
scheme presently under development. 
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Predicting the Structure of the Solar Corona During the 
11 August 1999 Total Solar Eclipse 


Zoran Mikic, Jon A. Linker, Pete Riley^ and Roberto Lionello 

Science Applications International Corporation, 10260 Campus Point 
Drive, San Diego, California 92121. U.S.A. 

Abstract. We describe the application of a three-dimensional magne- 
tohydrodynamic (MHD) model to the prediction of the structure of the 
corona during the total solar eclipse of 11 August 1999. The calculation 
uses the observed photospheric radial magnetic field as a boundary condi- 
tion. This model makes it possible to determine the large-scale structure 
of the magnetic field in the corona, as well as the distribution of the solar 
wind velocity, plasma density, and temperature. The density was used to 
predict the plane-of-sky polarization brightness prior to the eclipse. The 
prediction is compared with an eclipse image taken in Turkey. 


1. Introduction 

Total solar eclipses offer an excellent opportunity to observe coronal streamers. 
During a total solar eclipse the moon blocks the bright light from the solar 
disk, so that the faint light scattered by the solar corona, which is more than 
a million times fainter than the photosphere, becomes visible. During totality 
the structures that characterize the white-light corona become visible, including 
prominences, helmet streamers, polar plumes, and coronal holes. Observers who 
witness a total solar eclipse invariably report that it is a beautiful sight to behold. 

On 11 August 1999 a total solar eclipse occurred in Europe, the Middle 
East, and India. A partial eclipse was visible from the North-Eastern United 
States, and many parts of Europe, North Africa, and West Asia. (For an archive 
of past and future eclipse paths, see Fred Espenak’s NASA eclipse home page. 1 ) 
Fortunately, the viewing conditions were excellent in Harput, in Eastern Turkey, 
where one of us (ZM) successfully observed the eclipse. 

In this paper we describe the a method of predicting the structure of the 
solar corona using a magnetohydrodynamic (MHD) model. By using this model 
it is possible to calculate the density in the corona, and thereby to deduce 
the brightness of coronal structures, given the observed magnetic field in the 
photosphere (a quantity that is measured routinely by several ground-based 
observatories). A description of the MHD model, as well as its application to 
the prediction of the structure of the corona during eclipses, is given by Mikic 
et al. (1999). 


1 See: http : //sunearth . gsf c . nas a . gov/eclipse/eclipse . html 
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Previous comparisons of our results with coronal and interplanetary obser- 
vations have shown that the large-scale structure of the solar corona is largely 
determined by the photospheric magnetic field distribution. It is remarkable 
that an MHD model, which incorporates the .balance of magnetic, plasma, and 
gravity forces in equilibrium, can be used to estimate the large-scale distribution 
of plasma, magnetic field, and solar wind in the corona, using only the measured 
magnetic field in the photosphere. 

This model can be used to make a prediction of the structure of the corona 
prior to eclipse day by using synoptic magnetic field measurements taken from 
the preceding solar rotation. After the calculation is performed, the white-light 
polarized brightness of the corona can be simulated from the MHD solution 
by integrating the electron density along the line of sight in the plane of the 
sky. This image can then be compared with coronal images taken during the 
eclipse. We published our prediction on the World Wide Web on 4 August 
1999 (http://haven.saic.com/corona/modeling.html). Our prediction was 
updated and finalized on G August 1999, five days prior to the eclipse. 


2. The MHD Model 


Over the past decade we have developed a 3D MHD model of the corona and 
inner heliosphere. At present, we primarily use a polytropic MHD model for our 
large-scale calculations. This model, briefly described here, relies on a simple 
form of the energy equation (an adiabatic fluid with a reduced poly tropic index). 
A more sophisticated model incorporating a more realistic description of the 
energy flow in the corona (including thermal conduction parallel to the magnetic 
field, radiation loss, coronal heating, and Alfven wave acceleration) is under 
development (Mikic et al. 1999), and will be used in the future to improve our 
description of the corona and solar wind. 

In the polytropic MHD model, the coronal plasma is described by the fol- 
lowing equations: 


47T 

V x B = — J, 
c 

_ _ 19B 

VxE = ”ai-’ 
i __ T 

E + -v x B = 77J, 
c 

^ + V-(pv) = 0, 


p (^ + v • Vv) = x B - Vp + pg + V- («/pVv) , 


( 1 ) 

( 2 ) 

(3) 

(4) 

(5) 


^ + V- (pv) = -(7 - l)pV • v, (6) 

where B is the magnetic field, J is the electric current density, E is the electric 
field, p, v, p, and T are the plasma mass density, velocity, pressure, and tem- 
perature, g is the gravitational acceleration, 7 is the polytropic index, 77 is the 
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resistivity, and 1 / is the kinematic viscosity. The plasma pressure is p = 2nkT, 
where n is the electron and proton density, which are assumed to be equal. In 
the polytropic model, the complicated energy flow in the corona is modeled with 
a simple adiabatic energy equation with reduced 7 (i.e., smaller than 5/3; see 
Parker 1963). The primary motivation for using a reduced 7 is the fact that 
the temperature in the corona does not vary, .substantially (the limit 7 1 

corresponds to an isothermal plasma). A typical choice, used here, is 7 = 1.05. 

We have developed a three-dimensional code to solve equations (l)-(6) in 
spherical coordinates (r,0, <j>) (Mikic & Linker 1994; Lionello, Mikic, & Schnack 
1998). This code has been used extensively to model the 2D and 3D corona, 
including the structure of helmet streamers (Mikic & Linker 1996; Mikic, Linker 
& Schnack 1996; Linker et al. 1999), coronal mass ejections (Mikic & Linker 
1994; Linker, Mikic, & Schnack 1994; Linker & Mikic 1995; Mikic & Linker 
1997; Linker & Mikic 1997), and the long-term evolution of the solar corona 
and heliospheric current sheet (Mikic et al. 1999). Related techniques have 
been developed by Pneuman and Kopp (1971), Endler (1971), Steinolfson et 
al. (1982), Washimi, Yoshino, and Ogino (1987), Wang et al. (1993), Usmanov 
(1993, 1999), and Suess et al. (1999). 

A boundary condition is imposed 011 the radial component of the magnetic 
field, B r , at the lower corona (r = R s ). We use synoptic maps of the line-of-sight 
photospheric magnetic field measured at the National Solar Observatory at Kitt 
Peak (NSOKP) for this purpose. These photospheric magnetic field maps are 
built up from daily observations of the Sun during a solar rotation, and give a 
good approximation of the Sun’s magnetic flux if the large-scale flux does not 
change considerably over a rotation. We also specify a uniform plasma density 
(no = 2 x 10 8 cm“*) and temperature (To = 1.8 x 10 6 K) at r — R s in regions 
where the radial velocity is positive. 

The application of the model is described by Mikic et al. (1999). Briefly, 
we start with a potential field in the corona (with V x B = 0) that matches 
the measured field B r at the base of the corona, and a transonic spherically 
symmetric solar wind solution (Parker 1963) to specify p, p , and v. This initial 
nonequilibrium field is integrated in time until a steady state is reached. The 
final state has closed magnetic field regions (helmet streamers), where the solar 
wind plasma is trapped, surrounded by open fields (coronal holes), where the 
solar wind flows freely along magnetic field lines, accelerating to supersonic 
speeds. 


3. The Prediction 

Previous applications of our model to the prediction of the state of the solar 
corona were performed for eclipses that occurred close to solar minimum: the 3 
November 1994 eclipse 2 , which occurred during the declining phase of the last 
solax cycle; the 24 October 1995 and 9 March 1997 eclipses, which occurred at 
solar minimum, and the 26 February 1998 eclipse, which occurred during the 


2 We did not make a prediction for the 1994 eclipse. However, we compared eclipse images with 
a coronal model after the eclipse. 
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early rising phase of the new solar cycle. The 11 August 1999 eclipse presented a 
new challenge to our modeling efforts because" it occurred during the late rising 
phase of the solar cycle on the approach to solar maximum. During this time, the 
solar magnetic field is significantly more complex than during solar minimum. 
The number of active regions is much larger, and the streamer structure in the 
corona is considerably more complex. A consequence for numerical modeling is 
the need for a significantly larger number of mesh points in the calculation to 
resolve coronal structures. Whereas our solar minimum calculations were per- 
formed on (r, 6 , cj>) grids with 81 x 81 x 64 mesh points, in the present calculation 
we use 111 x 141 x 128 mesh points. In addition, the large-scale magnetic field 
strength is larger, requiring smaller time steps. These calculations therefore took 
significantly longer to perform (90 CPU hours on a single processor of the T90 
Cray supercomputer). We are presently developing a massively parallel version 
of our code which will enable us to perform these calculations in less time. 

Using magnetic field measurements from a previous solar rotation for a 
prediction is expected to be a poorer approximation during this phase of the 
solar cycle than during solar minimum, when the large-scale structure of the 
Sun changes slowly between solar rotations. The photospheric magnetic field 
evolves more rapidly, making synoptic magnetic field measurements a less reli- 
able approximation to the true state of the photospheric magnetic field. Figure 
1 shows the evolution of the NSOKP photospheric magnetic field synoptic maps 
during the four solar rotations preceding the eclipse. Note that there is signifi- 
cant evolution of the magnetic field from rotation to rotation as active regions 
emerge and disperse. Figure 2 shows a comparison of a NSOKP synoptic map 
from June-July, 1999 with one near solar minimum (August-September, 1996), 
showing that the photospheric magnetic field is considerably more complex now. 
Consequently, the coronal magnetic field would be expected to be significantly 
more complex also. This expectation is confirmed by our results, described 
below. 

The magnetic field is fitted at the solar poles using a smoothing proce- 
dure, since the line-of-sight component is not measured accurately there due 
to projection effects. Additionally, the data is smoothed everywhere, since the 
spatial resolution of the calculation is less than that of the measurements. The 
poor accuracy of the polar field (a limitation of the measurements) is expected 
to compromise our prediction more severely at this time of high solar activity, 
when active regions and coronal streamers are found at high latitudes, than at 
solar minimum. 

On 28 July 1999 we started an MHD computation to predict the struc- 
ture of the solar corona during the 11 August 1999 eclipse. We used NSOKP 
photospheric magnetic field measurements from Carrington rotation 1951 (corre- 
sponding to the dates 24 June-21 July, 1999). Figure 3 shows the raw NSOKP 
data and the smoothed data used in our calculation. It is apparent that the 
large-scale features in the magnetic field are reasonably well represented by the 
smoothed field. 

Once the coronal solution reaches equilibrium, we use the density from the 
calculation to simulate the white-light polarized brightness of the corona on 
eclipse day. The polarized brightness is computed by integrating the electron 
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Evolution of the Photospheric Magnetic Field 
Prior to the 11 August 1999 Total Solar Eclipse 

Carrington Rotation 1949 (Mav 1 - May 28, 1999) 




Figure 1. Synoptic maps of the radial component of the photospheric 
magnetic field for three solar rotations (Carrington rotations 1949- 
1951), as measured at the National Solar Observatory at Kitt Peak. 
Black shows fields directed into the Sun, whereas white shows fields 
directed out of the Sun. These maps show that the magnetic field is 
changing rapidly during the approach to the eclipse of 11 August 1999. 
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Comparison of Photospheric Magnetic Fields: 

Near Solar Minimum and Approaching Solar Maximum 



Figure 2. Synoptic maps of the radial component of the photospheric 

magnetic field near solar minimum (top) and approaching solar maxi- 
mum (bottom). 


density along the line of sight in the plane of the sky (convolved with a scattering 
function and filtered with a radial “vignetting function”). 

On 6 August 1999 we performed an updated calculation using the most 
recent magnetic field measurements available from NSOKP on that day, to in- 
corporate the evolution of the magnetic field. The measurements used in the 
synoptic map were those from 13 July-5 August, 1999 (the synoptic map in- 
cluded Carrington longitudes 0°-107° from Carrington rotation 1951 and 107°- 
360° from rotation 1952). Due to time constraints, the updated calculation was 
performed on a coarser grid (61 x 71 x 64), and required 4 hours of CPU time. 
The results show that there are changes in the predicted coronal structure in 
the updated calculation. In particular, the position of some streamers changed 
perceptibly, and they* had different tilts (especially near the west solar limb, 
where the magnetic field data was updated m6st significantly). However, the 
large-scale structure of the corona did not change dramatically. The updated 
calculation was not relaxed to steady state sufficiently, and the spatial resolution 
was inadequate, so that the results shown in this paper use the calculation with 
data for Carrington rotation 1951. 
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Photospheric Magnetic Field for Carrington 
Rotation 1951 (June 24 *■ July 21, 1999) 

Kitt Peak SvnoDtic Mao 



Figure 3. Synoptic map of the radial component of the photospheric 
magnetic field measured at the National Solar Observatory at Kitt Peak 
for Carrington rotation 1951 (top), and the smoothed version (bottom) 
that was used as a boundary condition for the MHD calculation for the 
eclipse prediction. 


4. Comparison with an Eclipse Image 

NASA astronomer Fred Espenak photographed the total solar eclipse on 11 
August 1999 from Lake Hazar, Turkey. Twenty-two separate images with differ- 
ent exposures were combined digitally into a single composite image that more 
closely resembles the appearance of the solar corona as seen by the human eye 
(Espenak 1999). Additional photographs of the eclipse can be found on Fred 
Espenak’s eclipse World Wide Web site. 3 

Figure 4 shows this eclipse image, as well as our predicted polarization 
brightness (pB ) from the MHD model. The pB image was calculated from the 
simulated corona as it would appear in the plane of the sky on 11 August 1999, 
at 11:38 UT (corresponding to totality in Eastern Turkey), and has been radially 
detrended to account for the fall-off of coronal brightness with distance from the 


3 See: http : //vww . mreclipse . com/TSE99report s /TSE99Espenak . html 
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Comparison of a 3D MHD Coronal Prediction with 
an Image of the 11 August 1999 Total Solar Eclipse 



Predicted Polarization Brightness (MHD Model) 



Predicted Magnetic Field Lines (MHD Model) 



Figure 4. Comparison between a composite eclipse image created 
from photographs taken by Fred Espenak in Lake Hazar, Turkey (top) 
with the predicted polarization brightness of the simulated solar corona 
from our 3D MHD model (middle). The projected magnetic field lines 
from the model are also shown (bottom). Terrestrial (geocentric) north 
is vertically upward. (Solar north is 15° counterclockwise.) The eclipse 
image is copyrighted © 1999 by Fred Espenak. 
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Sun. This detrending is similar to the effect of combining various exposures 
in the composite eclipse image (although it is not an identical procedure), and 
is an attempt to reproduce the image seen by the human eye. Therefore, the 
agreement between- the two images can only be expected to be qualitative at 
best; namely, the position and structure of streamers and coronal features ought 
to agree, but the brightness of individual features would not be expected to 
agree. 

It is apparent that the eclipse image has considerably more fine structure 
than our prediction. This is primarily due to the limited spatial resolution of 
our calculation, but may also be due to the polytropic energy equation used 
in our model, and the fact that the fine-scale structure has been emphasized 
in the composite image (Espenak 1999). Note that the features that resemble 
“plumes,” that are visible in the eclipse image at equatorial latitudes on the east 
solar limb, are not present in the model. 

Previous comparisons of our simulations with Mauna Loa MK3 corona- 
graph observations on several days surrounding our solar minimum eclipse pre- 
dictions have confirmed that the basic large-scale three-dimensional structure 
of the streamer belt is captured in our model (Linker et al. 1999). However, 
the agreement between our prediction and eclipse images as we approach solar 
maximum is not as good. The principal discrepancy seems to be in the tilt of the 
helmet streamers, and in the fact that we miss individual streamers, especially 
in the polar regions. Possible reasons for this discrepancy are discussed in the 
next section. Overall, however, the prediction gives a fair approximation to the 
state of the large-scale solar corona. 

It is difficult to portray the complex three-dimensional structure of the 
solar corona for this simulation using still images. It is best seen in a movie of 
the polarization brightness during the whole solar rotation (Carrington rotation 
1951), which can be found at our Web site. 4 


5. Discussion 

During solar maximum the solar magnetic field is changing rapidly as new active 
regions emerge daily and disperse as they interact with existing active regions. 
The accuracy of an equilibrium model, which is based on synoptic magnetic field 
measurements, such as the one we describe, is therefore expected to be limited. 
Not surprisingly, the match between the prediction and actual observations is 
poorer during this time of high solar activity than it was during solar minimum. 
Furthermore, the failure to match the transverse component of the magnetic 
field (which is presently not measured), which specifies the “shear” in the field, 
and consequently the energization level of the coronal field (Mikic & Linker 
1997), is expected to be most severe at solar maximum when the magnetic field 
is generally most active (and hence significantly “energized”). 

These two key approximations will be improved in future versions of our 
model. First, we now have the capability of running time-dependent simula- 
tions of the evolving solar corona as it responds to changes in the photospheric 


4 See: http : //haven . saic . com/ corona/modeling . html 
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magnetic field. In the future, we will be able to use daily magnetic field mea- 
surements to calculate the evolving solar corona.' Second, current plans at Kitt 
Peak National Solar Observatory call for measurements of the vector magnetic 
field in the photosphere to be taken (the SOLIS project), which can be used 
in our model to match the transverse component of the magnetic field at the 
base of the corona (Mikic et al. 1999). A benefit of vector magnetic field mea- 
surements will be that polar fields will be more accurately determined (since 
the polar field, which is predominantly radial, will be measured as part of the 
transverse component). This can be expected to improve future predictions. 
Finally, our predictions will also be improved when our model with the more 
sophisticated energy equation becomes operational, and when we can perform 
higher-resolution calculations in less time with the massively parallel version of 
our code. 
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Abstract 

The magnetic field dominates the structure of 
the solar corona. We have developed a magne- 
tohydrodynamic model of the large scale struc- 
ture and dynamics of the corona with a self- 
consistent treatment of the thermodynamics 
that includes radiation losses, thermal con- 
duction along the magnetic field, and coro- 
nal heating. We specify a magnetic flux dis- 
tribution on the solar surface and integrate 
the time dependent MHD equations to steady 
state. The model has been used to simulate 
the streamer/coronal hole configuration in two 
dimensions from the upper chromosphere and 
transition region to 30 fZ©. 

1 Introduction 

The large scale structure and dynamics of the 
solar corona is dominated by the magnetic 
field. In the past the hydrodynamics and ther- 
modynamics of the closed field regions (loops) 
of the corona has been studied with ID models 
by many, amongst whom are [9], [5], [12], [10]. 
Other ID models such as [14] and [3]. have 
been used to reproduce the solar wind in the 


open field regions. [8] successfully developed 
a 2D magnetohydrodynamic (MHD), isother- 
mal model of the corona with a streamer and 
the open field regions. A 3D MHD, poly- 
tropic model could reproduce many features of 
the solar corona during Whole Sun Month [6]. 
However, since a simplified energy equation is 
used, it is not surprising that those models 
do not show the contrast in the plasma den- 
sity and temperature between the open and 
closed field region and the contrast in speed 
between the fast and slow solar wind. [11] 
and [13] used a 2D MHD model of the corona 
and of the solar wind which includes heat- 
ing and thermal conduction. In this paper 
we present a self-consistent study of the 2D 
large scale structure and dynamics of the so- 
lar corona and of the upper chromosphere. It 
Jias been obtained with an MHD numerical 
algorithm which treats the thermodynamics, 
including thermal conduction parallel to the 
magnetic^ field, radiation, and coronal heat- 
.ing. Thermal conduction is collisionally dom- 
inated in the inner corona, where it is calcu- 
lated with Spitzer’s formula, and smoothly be- 
comes collisionless in the outer corona. The 
results show sharp gradients in the transition 
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region, which appears structured according to 
the field topology. The temperature and den- 
sity contrasts between the streamer and the 
poles are evident too. The aim of the present 
investigation has not been to reproduce the 
speed contrast between the fast and slow wind, 
which will be addressed with the same model 
at a later time. 

2 The MHD and Thermo- 
dynamic Model 

We have used our computational model to sim- 
ulate the large scale structure of the 2D solar 
corona and the transition region in spherical 
coordinates (r, 6). Our model is capable to 
perform 3D simulations as well. It includes 
a self-consistent treatment of the magnetic, 
thermal and gravitational forces and of the en- 
ergy transport. The physics is described in 
details in [7]. We solve the resistive and vis- 
cous MHD equations in spherical coordinates. 
The energy equation includes a radiation law 
function as in [1], thermal conduction along 
the magnetic field and coronal heating. Ther- 
mal conduction is collisional (Spitzers law) for 
R© < r < 10 R© and smoothly becomes colli- 
sionless [4] for r > 10 R©. 

3 The Simulation 

We have started our simulations by imposing 
an initial condition extracted from a ID so- 
lution of a solar wind which includes the up- 
per chromosphere and the transition region. A 
dipole field of 2.205 Gauss at the poles is pre- 
scribed as initial magnetic field. At the base 
of the domain, r = 1 R©, the magnetic flux 
corresponding to the aforementioned dipole is 
held fixed during the simulation. A fixed tem- 
perature of 2 x 10 4 K and a fixed density of 
1.5 x 10 10 cm" 3 are also prescribed at the in- 
ner boundary. These values are compatible 
with those of the upper chromosphere* Char- 


acteristics equations are used to calculate the 
velocity along the magnetic field. These are 
also employed at the outer boundary (30 R©), 
which is placed beyond the sonic and Alfvenic 
points. Because of the uncertainties surround- 
ing coronal heating mechanisms, the heating 
in the energy equation is chosen to be a pa- 
rameterized function 

H ch = H 0 exp[-(r-Ro)/X\. (1) 

In the present simulation heating has been 
specified uniform in 8 and A = 0.7 R©. Hq 
is such that the heat flux at the base is 10 5 
erg cm -2 s” 1 . The Lundquist number we have 
used in the calculation is 5 = 5 x 10 4 , which 
is many order of magnitude smaller than the 
value in the corona. This is an inevitable 
drawback common to all numerical compu- 
tations, which are limited by the grid res- 
olution. The same argument applies in or- 
der to explain the choice of a finite viscos- 
ity z/, which corresponds to a viscous diffusion 
time r u = 3.33 x 10 3 ta , where r v = R^jv, 
ta — Rq/ v a is the Alfven transit time, and va 
is the Alfven speed. We have used a 201 x 301 
nonuniform grid in (r, 9). Radial grid points 
are denser at the base of the domain, the ra- 
tio between mesh spacing at 30 R© and 1 R© 
being 2 x 10 4 . The ratio between latitudinal 
mesh spacing at the pole and at the equator is 
2. We have advanced the MHD equations and 
looked for a steady solution. 

4 Results 

During the simulation a considerable fraction 
of the closed magnetic lines of the initial dipole 
are progressively opened. The final configura- 
tion of the magnetic field features a helmet 
streamer with a cusp-like top surrounded by 
open field regions (see Fig. 1). A current sheet 
is formed on top of the streamer, as Fig. 2 
shows. 

Currents are present in the closed field re- 
gion under the streamer as well. This is in 
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Figure 1: Plasma temperature at increasingly smaller scales with superimposed magnetic field 
-lines. In the left corner, the temperature scale. In the right corner a Cartesian projection. 
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Figure 2: Current density contours at the end 
of the simulation. The helmet streamer and 
the heliospheric current sheet are visible. 

contrast with the results obtained by [8], who 
found currents only at the boundary separat- 
ing the streamer from the open field regions. 

In Fig. 3a we plot the radial velocity as a 
function of the radius for the pole and the 
equator. The contrast between the two is in- 
different. We have not included in this simula- 
tion a heating term dependent on the latitude 
and an additional source in the momentum 
equation, which are necessary to reproduce the 
fast and slow wind [13]. In fact the aim of this 
investigation is to study the 2D structure of 
the upper chromosphere, the transition region 
and the corona, and not to reproduce all the 
characteristics of the solar wind. Figures 3b 
and 3c show the particle number density and 
the temperature. A jump of several order of 
magnitudes is evident in the transition region 
in both cases. Notice that the equatorial re- 
gion has higher temperature and density in the 
radial range corresponding to the streamer. 

In Fig. 4a we plot three magnetic field lines 
representing a low lying loop, an intermediate 
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Figure 3: Radial velocity (a), particle num- 
ber density (b), and temperature (c), as func- 
tions of the radial distance in solar radii for 
the equatorial and polar stream. 
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Figure 4: Three field lines (a) representing a 
low lying loop, an intermediate one, and the 
longest loop present in the simulation (on the 
solar surface the magnetic flux distribution is 
showed). Number densities (b) and tempera- 
ture (c) as functions of the distance measured 
along the field lines. 



Figure 5: Polarized brightness image obtained 
from the data of the simulation. 


loop, and a loop at the boundary between the 
open field region and the streamer. On the 
solar surface we show the magnetic flux dis- 
tribution. Figure 4b has the densities for the 
three loops as function of the distance mea- 
sured along the field lines. The temperature 
for the same loops is presented in Fig. 4c. 

Figure 1 shows the temperature at increas- 
ingly smaller scales with superimposed mag- 
netic field lines. The temperature maximum 
of about 2 x 10 6 K is in the closed field re- 
gions under the streamer, where hot plasma is 
trapped by the magnetic field. A projection in 
Cartesian coordinates shows that the height at 
which the transition region is found appears to 
be dependent on the latitude. 

From our data we have produced a syn- 
thetic polarized brightness image (Fig. 5). The 
features can be compared with those present 
in real observations, taking into account that 
there is no occulting disk in our image. 
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5 Conclusion 

We have been able to perform a two dimen- 
sional MHD simulation of the solar corona, 
the transition region and the upper chromo- 
sphere. This has been possible since our code 
includes a self consistent treatment of the ther- 
modynamics with parallel thermal conduction, 
heating and radiation losses. It has been a 
very challenging computation, as the temper- 
ature, density and velocity span several orders 
of magnitude across very small scales (~ 10 -3 
iZ 0 ). The simulation has allowed us to inves- 
tigate the large scale structure of the coro- 
nal magnetic field and underline the connec- 
tions with the plasma properties. We have 
reproduced the helmet streamer surmounted 
by the heliospheric current sheet. The closed 
field regions below the helmet streamer appear 
to contrast in temperature and density with 
the open field regions, through which the so- 
lar wind flows. This seems to justify the pres- 
ence of the currents in the closed field regions, 
where a Lorentz’s force balancing the pressure 
gradients is needed. This is in contrast with 

[8], who however used the isothermal approx- 
imation in their model. 

The transition region appears to be struc- 
tured differently according to the magnetic 
field topology. In particular the height at 
which the temperature gradient rises is lower 
for the longest and hottest loops in the 
streamer than for the open field region and 
the shorter loops. This structure is also visible 
in a polarized light image, where the streamer 
gives the impression to possess two "legs" . 

With the present model we plan to study the 
2D and 3D structure of the corona and of the 
solar wind. The latter can be modeled by pre- 
scribing a nonuniform heating functions and 
an additional momentum source in the calcu- 
lation. Since our model includes the thermo- 
dynamics now, it may be possible to employ it 
to simulate prominence formation. 
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